Exact Study of the ID Boson Hubbard Model with a Superlattice Potential 



O 
O 

(N 

X 



5-1 

43 



V.G. Rousseau, 1 D.P. Arovas, 2 M. Rigol, 1 F. Hebert, 3 G.G. Batrouni, 3 and R.T. Scalettar 1 

'Physics Department, University of California, Davis, California 95616, USA 
2 Physics Department, University of California, San Diego, California 92093, USA 
Institut Non Lineaire de Nice, 1361 route des Lucioles, 06560 Valbonne, France 

We use Quantum Monte Carlo simulations and exact diagonalization to explore the phase diagram 
of the Bose-Hubbard model with an additional superlattice potential. We first analyze the properties 
of superfluid and insulating phases present in the hard-core limit where an exact analytic treatment 
is possible via the Jordan- Wigner transformation. The extension to finite on-site interaction is 
achieved by means of quantum Monte Carlo simulations. We determine insulator/superfluid phase 
diagrams as functions of the on-site repulsive interaction, superlattice potential strength, and filling, 
finding that insulators with fractional occupation numbers, which are present in the hard-core case, 
extend deep into the soft-core region. Furthermore, at integer fillings, we find that the competition 
between the on-site repulsion and the superlattice potential can produce a phase transition between 
a Mott insulator and a charge density wave insulator, with an intermediate superfluid phase. Our 
results are relevant to the behavior of ultracold atoms in optical superlattices which are beginning 
to be studied experimentally. 

PACS numbers: 05.30.Jp,71.10.Fd,02.70.Uu 
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I. INTRODUCTION 

Originally proposed as a model for short-coherence 
length superconductors and Josephson junction arrays, 
the boson Hubbard Hamiltoniar. 1 has over the last sev- 
eral years been widely used to understand the physics of 
ultra-cold optically-trapped atoms; 2 While previous work 
had considered translationally invariant systems 3 ^«yi*L& 
or the effect of a random chemical potential in order 
to understand glassy behavior a key feature of the 
model as applied to cold atoms is the inclusion of a (usu- 
ally quadratic) external potential which reflects the mag- 
netic confinement. This potential leads to a number of 
interesting effects including the coexistence of superfluid 
and Mott insulating regions within the trapi2*iiii2i 

Early on in the application of the Bose-Hubbard model 
to optically confined atoms, the study of a "superlattice" 
in which the confining potential has multiple minima was 
considered^ Upon increasing the chemical potential, the 
density profile evolves from a situation where the boson 
densities in the different minima are independent to one 
where superfluid 'necks' develop and join the bosons in 
the different minima. Subsequently, further mean field 
theory treatments developed a more quantitative under- 
standing of the phase diagram, in the case when the su- 
perlattice potential varies with a period of T = 2,3 and 
4 sites lii^i^ii^ Mott insulating phases with fractional fill- 
ings exist, and, interestingly, under certain circumstances 
the usual Mott phase at p = 1 can be absent. The 
physics of the Bose-Hubbard Hamiltonian with aperiodic 
potentialsjiSiil in which localization without disorder can 
occur, has also been discussed^ 

Various experimental realizations of such multiple well 
superlattices have been proposed, from a double well 
magnetic trap for Bose-Einstein condensations in which 
the barrier height and well separation are smoothly 
controllable, 19 to periodic potentials where the lattice 



constant is especially large, allowing the loading of many 
bosons per minimum^SiSi Superlattice potentials similar 
to the ones considered in this work have been realized by 
Peil et al22k One of the great advantages of these ultra- 
cold gas realizations of strongly correlated systems is the 
experimental possibility to tune all parameters at will. 

In this paper we use Quantum Monte Carlo (QMC) 
simulations and exact diagonalization to study the 
physics of the boson Hubbard Hamiltonian and its in- 
finite U limit in the presence of a superlattice potential. 
In contrast to previous mean-field studies ji&iaiili we con- 
sider the hopping parameter to be independent of the 
position in the lattice. In addition, our numerical ap- 
proaches provide an exact treatment of correlations that 
are particularly important in one dimension. We also 
compute important quantities like the superfluid density 
and the momentum distribution function, which make 
connections with experiments. Superlattice potentials 
similar to the ones considered in this work have been re- 
alized by Peil et alml and Sebby-Strabley et a/.— One of 
our main results is that the superlattice produces insulat- 
ing phases for commensurate fractional fillings. For the 
hard-core case, this behavior can be explained in terms 
of band structures by performing an exact mapping onto 
a spinless fermionic system. Insulating behavior persists 
for the soft-core case with sufficient on-site repulsion. On 
the other hand, at integer fillings, changing the ratio be- 
tween the on-site repulsion and the strength of the ad- 
ditional superlattice potential can produce an intermedi- 
ate superfluid phase between Mott insulating and charge 
density wave phases. 

The exposition is organized as follows. In Sec. |H] we 
study the infinite U limit of the boson Hubbard Hamil- 
tonian with an additional superlattice potential. The 
generalization to finite on-site interactions is presented 
in Sec. 11111 where we use QMC simulations to solve the 
problem exactly. We also present in Sec. IIIII an analy- 
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sis of the atomic limit that helps to understand results 
for large but finite U, where multiple occupancy of the 
lattice sites is allowed. In Sec. IIVI we present a theo- 
retical analysis extending the atomic limit by allowing 
finite hopping amplitudes, which gives further insight on 
our numerical results for large values of U . Finally, the 
conclusions are presented in Sec. [3 

The translationally invariant boson Hubbard modeli 

is: 

H = -t^{a\ aj + ojoi) + U^hi{fn - 1) (1) 

<«> 

The operators a\,ai create (destroy) a boson on site i, 
and obey commutation rules [ffli, at] = 6ij. The number 

operator is hi — a\di. The hopping parameter t measures 
the kinetic energy and U the strength of the on-site re- 
pulsion. We will consider a one dimensional lattice so 
that the sum (ij) over near neighbors has j = i + 1. The 
ground state phase diagram is well known pL&iSiSiSi At 
commensurate fillings, and for sufficiently large U, the 
bosons are in a gapped Mott insulating phase. Away 
from integer filling, or for weak coupling, the system is 
superfluid. In the limit U — ► oo this model maps onto 
the spin 1/2 XY model, with the z component of mag- 
netization playing the role of the boson density. 

To obtain a superlattice, we consider a case where a 
low amplitude, long wavelength potential is added to the 
usual high intensity short wavelength optical potential 
which generates the lattice in which the atoms move. 
(See Fig. Atoms in the resulting superlattice thus 
have a hopping parameter which is independent of spa- 
tial position. This is completely analogous to the usual 
optical trap configuration and associated model calcu- 
lations. The long period potential we consider has the 
form 

Kxt = A^cos^ hj. (2) 

3 

We will be interested in understanding the ground state 
phase diagram as a function of the energy scales U/t, A/t, 
particle density p, and the period T. 



II. ANALYTIC TREATMENT OF THE 
HARD-CORE LIMIT 

We first consider the hard core limit, U — oo, which is 
exactly solvable via the Jordan- Wigner transformation, 

«j = c ] n ^ . 0) 

n=l 

If the boson operators obey general twisted boundary 
conditions a f L+1 = e~ lS a\, it proves convenient to invoke 
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FIG. 1: (Color online) An example of a superlattice potential 
we study in this work. The long period potential provides 
a small additional modulation to the deep optical potential 
which produces the lattice. As a consequence, the hopping 
parameter t is independant of position. 

a unitary transformation to transfer the boundary con- 
dition phase from the bosonic operators to the hopping 
integrals, viz. 

L L f2im\ 

H = (Cn+ia« a n+l+ H - c -) +yl X! C0S ( ~ ) a n a n, 

n=\ n=l ^ ' 

(4) 

where = a\ and t' n n+1 = t for 1 < n < L and 

t' L 1 = e iS t. When mapping this Hamiltonian, with the 
additional on-site hard-core constrains 

4 2 = a « = 0, {a n ,at} = l, (5) 

onto a non-interacting fermionic Hamiltonian one notices 
that 

L 

a\a L = -c\c L ft e- c K (6) 

n=l 

which means that the equivalent fermionic Hamiltonian 
takes the form 

L L /27rn\ 

H = -^2 (t n< n+lcicn+i+R.C.)+A^2cosl — Jc+Cn , 
n=l n=l ^ ' 

(7) 

where c\ +l — c\ and t n ,n+i = t for 1 < n < L and t^ i — 

e i{&+ri) ^ depending on whether N = (5^n=i c n Cn ) IS even 
(rj = 7r) or odd (77 = 0). Via a second unitary transforma- 
tion, we can impose the boundary phase uniformly over 
the hopping integrals, yielding t n ,n+\ — e^ s+v ^ L t for all 
n. 

We now transform to a quasi- momentum basis, writ- 
ing cjj = L~ x l 2 e~ lkn c' k , where k is quantized with 
e %kL _ ^ xhe superlattice potential couples states of 
quasi-momenta k and k±Q, where Q = 2tt/T. Restrict- 
ing k to the reduced Brillouin zone (BZ) [ — the 
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FIG. 2: (Color online) Band structure of the hard-core Bose 
Hubbard model for periods T — 2,3, 4, 6. The two central 
energy bands in the T — 4 case touch, and the system is not 
an insulator at p = | . 




FIG. 3: (Color online) The gap A as a function of T for 
t = 1 and A = 2 at half filling. A is finite only for T = 4p + 2. 
The inset shows the same plot but with logarithmic axes, 
emphasizing a power law decay which is discussed in the text. 



matrix elements of TL are 

«,*,, = (k + lQ\H\k + l'Q) 

= -2*cos(fc + C + lQ)$' + ^Af , (8) 

where £ = (5 + il)/L, and (5 ; T ; , = <5u'modT- This defines 
a T x T matrix for each A; in the reduced zone. The 
eigenvalues give the T energy bands E a (k). The T = 2 
case is familiar: 



E±(k) = ±y/4t 2 cos 2 (k + ()+A 2 



(9) 



The effect of the periodic potential is to open up gaps 
at the boundaries of the reduced BZ. Figure [2] shows the 
band structure for T = 2,3, 4, 6. There are energy gaps at 
i for T = 2, at ± and § for T = 3, at p = \, §, §, §, § 
6. However, for T = 4 the energy bands cross at 
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p = i and gaps exist only at i and | . 

Taking t/A — * yields the atomic limit. For t = 0, the 
energy levels are = A cos(2%n/T). For small i/A, the 
bandwidth of each of the T bands may be obtained to 
leading order by appealing to the locator expansion for 
the Green's function, as we show in Appendix A. 

Because the number of bands is equal to the number of 
sites per supercell, one expects that the system is insulat- 
ing at half filling when the period T is an even number. 
For such cases the first T/2 bands are completely filled, 
and the creation of an exciton requires a finite supply of 
energy. However it can happen that the valence band 
crosses the conduction band, as illustrated in Fig. for 
T = 4. Then no gap arises, and the system is not insu- 
lating. As shown in Appendix A, this situation occurs 
(for p = ^) whenever T — 4p, p being an integer. The 
situation is depicted in Fig. |21 For T = Ap + 2 a gap 
occurs, and decays as a power law as a function of T for 
a fixed value of A. (This is shown in the inset, where a 



fit provides an exponent of —1.03 for A/t = 2). In the 
atomic limit, we have A = 2Asm{^), hence A cx T _1 
at large T. For odd values of T, the system at half- filling 
is never insulating because the highest occupied band is 
itself half-filled. 

The dimensionless superfluid density is given by the 
expression 



L d 2 F 
2tW 
Tiv-p 

27rf 



5=0 



(e = o,l = 00) 



(10) 

(11) 



where F is the free energy and vp = j- ^\k F ^ s ^ e F errm 
velocity. The second expression, valid at temperature 
= in the thermodynamic limit L = 00, follows from 
the fact that 8E n /d5 = L^ 1 dE n /dk. We see from this 
last expression that the superfluid density vanishes when- 
ever the Fermi level lies within a band gap. Furthermore, 
if the last partially occupied fermion band is nearly filled 
or nearly empty, then vf ~ hk-p/m*, with the effective 

mass (m*)" 1 = ^j-^pS and the dimensionless density is 
p — pMott + &Pi with Sp = ±fejr/7r. Therefore we find 



Ps(p) 



jOMott 



(12) 



where m = h 2 /2ta 2 is the bare 'mass' (a being the physi- 
cal lattice constant), and pMott is the density correspond- 
ing to an integer number of filled bands. 

Another definition of superfluid density, which be- 
comes exact in the thermodynamic limit, is based on the 
free energy difference between periodic (6 — 0) and an- 
tiperiodic (<5 = n) boundary conditions: 
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(13) 
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20 lattice sites 
100 lattice sites 
Thermodynamic limit 



T=10 p=0.25 




FIG. 4: (Color online) The superfluid density, p s , with 
p = 0.25 and T — 10 as a function of A. We show results 
for two values of the system size L, and the thermodynamic 
limit result Eq. Ijlljl . No difference is perceptible between the 
results for L = 100 and L — oo. 
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FIG. 5: (Color online) p s L and gap A for T = 10 at half 
filling, as functions of A. For this filling, the critical value 
of A for which p$ vanishes shows a tendency to go to zero, 
when increasing the size of the system. This is due to the gap 
A which vanishes only at A — 0. For small but finite A, A 
follows a power law with the exponent 4.97. 



where Fg is the free energy at boundary phase 5. Eq. 
is of course equivalent to a discrete approximation to the 
second derivative in Eq. I|10[) . but it has the advantage 
of being easier to evaluate numerically. We find that 
in the superfluid phases finite size effects related to the 
definition in Eq. i|13|) start to be negligible for quite small 
system sizes, of the order of L w 100. This is shown in 
Fig. 0] which displays the values of the superfluid density, 
p s , for p = 0.25, T = 10 as a function of A obtained 
for two finite systems (L = 20, 100) following Eq. I|13|) 
and the thermodynamic limit calculation following Eq. 
pi[l. The curve for L = 100 is perfectly superposed to 
the thermodynamic limit result, and even L = 20 gives a 
very good estimation of p s . 



We should stress, however, that in the insulating 
phases, i.e., for commensurate fillings where p s is zero, 
finite size effects become relevant in Eq. ^ 1 of) when the 
gap approaches zero as A decreases. For those cases one 
needs to take the limit L — > oo to obtain the correct 
p s = 0. This can be seen in Fig. [S] where pf? [Eq. (|13|) ]. 
and the gap A, are plotted as functions of the amplitude 
A of the modulating potential, for T — 10, half filling, 
and different values of L. The inset shows the gap with 
logarithmic axes, emphasizing a decay with a power law 
when A goes to zero, with an exponent of 4.97, consistent 
with the exact value of 5 obtained from locator expansion 
of the Green's function for the Hamiltonian in Eq. JHJ). 

From the viewpoint of the fermion Hamiltonian, we 
might regard the insulating phases at fractional densities 
corresponding to complete filling of each E n (k) as 'band 
insulators'. However, when the hard core constraint is 
lifted, the insulating phases of the bosons are more prop- 
erly regarded as 'Mott insulators'. Consider, for example, 
the T = 2 case. The eigenvectors of the lowest energy 
band E_{k) have their largest density on those spatial 
sites with superlattice potential —A. If U is sufficiently 
weak, multiple occupancy of this spatial sublattice will 
not be energetically forbidden, and there is no reason to 
expect then an insulating phase at p = i. Thus, as U 
is decreased, one eventually will reach a quantum critical 
point where the gap vanishes and the system becomes su- 
perfluid. This is indeed what emerges from our Quantum 
Monte Carlo analysis of the soft-core model, discussed 
further below. 

In Fig.El we plot the superfluid density p s and energy 
gap A for the T = 10 system versus the filling p. As 
expected, A is finite only for p — j/T, where an integer 
j number of bands are completely filled. The superfluid 
density p s exhibits local maxima at p = (j + |)/T, in 
the centers of the bands. We note that p s is greatest in 
the central bands, and smallest in the outer bands. In- 
deed, the superfluid density must be small when p itself 




FIG. 6: (Color online) Superfluid density p a and gap A of 
the hard-core boson Hubbard Hamiltonian for T = 10, t = 1, 
and A — 2. Gapped band insulating phases exist at fillings 
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FIG. 7: (Color online) Superfluid density p s for T — 10 and 
different fillings, as functions of A, with t = 1. p s decays as a 
power law (see text for the exponents). 



is small, since there are fewer bosons to contribute to p s . 
As p increases (restricting our attention to the band cen- 
ters), p s increases concomitantly, until one passes p = |, 
when the trend reverses. For large p one should think in 
terms of superfluidity of holes (i.e. empty sites). 

Figure shows p s for several fillings as functions of 
A, for T = 10. The inset corresponds to the same plot 
but with logarithmic axes, showing a power law decay of 
p s . The numerically extracted exponents are —8.47 for 
p = 0.05, -6.94 for p = 0.15, -7.06 for p = 0.25, and 
—5.00 for p — 0.35 and p — 0.45. These values compare 
well with the exact results for to* derived earlier from the 
locator expansion: (to*) -1 cx t T jA T ~ x for bands arising 
from atomic levels = A cos (^f?) which are nondc- 
generate (n = 5 for p = 0.05), and (to*) -1 oc t 2n /A 2n ~ 1 
arising from degenerate atomic levels (n — 4 for p = 0.15 
and p = 0.25; n = 3 for p = 0.35 and p = 0.45). 

The full phase diagram for T — 10 is shown in Fig. |H1 
The gapped insulating phases extend all the way down 
to A = 0, as there is a finite gap for any finite value of 
A (Fig. ISJ) . As A/t increases, off the magic densities, p s 
goes to zero as a power law (Fig. 0), and for this reason 
the associated regions of the phase diagram (red color) 
are labeled 'weakly superfluid'. 

To conclude this section we discuss how the superfluid 
and Mott insulating phases introduced above could be 
detected in experiments with ultracold bosons on optical 
lattices. For that we study the imprint of these phases 
in the hard-core boson momentum distribution function 
[n(k)]r& which is a quantity that can be easily obtained 
in time of flight measurements. 

Figure shows n(k) for T = 10, p = 0.45 (superfluid 
case), and three values of A. For A = (top panel), a 
single peak can be seen at n(k = 0) with a height that 
scales proportional to the square root of the number N 
of particles in the system, and thus diverges in the ther- 
modynamic limits This peak signals quasi-long range 
one-particle correlations typical of the superfluid state in 
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FIG. 8: (Color online) Phase diagram of the hard-core Bose 
Hubbard model as a function of filling p and strength A of the 
superlattice potential. Fillings p = , , . . . are insulat- 
ing for all A 0. At fillings which are not commensurate with 
the superlattice potential, the superfluidity becomes small at 
large A. (See Figure |7J. 



ID. The introduction of the additional potential in Eq. 
(0) introduces a modulation in the one-particle correla- 
tions but does not destroy their quasi- long range order. 
As a consequence one can see in the two lowest panels of 
Fig. that additional sharp peaks appear with momenta 

U L "7T 

10 ' 

On the other hand, in the insulating phases (at comen- 
surate fillings) the opening of a gap in the one parti- 
cle excitation spectrum produces an exponential decay 
of the one-particle correlations. As seen in Fig. ^3 this 
exponential decay destroys the sharp peaks observed at 
n(k — 0) in the superfluid state (top panel). As a con- 
sequence, a very broad distribution is observed in n(k) 
with no additional satellite peaks at k = ±^ (two low- 
est panels). Hence, measuring n(k) in experiments would 
unambiguously differentiate between the superfluid and 
insulating phases we have analyzed in this section. 



III. QUANTUM MONTE CARLO 
SIMULATIONS OF THE SOFT-CORE CASE 

As a first step in studying the soft-core case, it is worth- 
while to recall the phase diagram of the uniform case 
A = in the (p/U, t/U) plane fFig.lTTl). In this situation 
the system is superfluid, except for integer densities and 
sufficiently large on-site repulsion. In this latter case, we 
have an incompressible Mott insulator. The transitions 
between superfluid and Mott insulating phases in the d- 
dimensional boson Hubbard model are known to be mean 
field like when driven by a change of the density, and of 
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Mott region with p = 
for p = 2 starts at U 
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FIG. 9: The momentum distribution function for p = 0.45. 
Top panel : The superlattice is turned off, and we 
have a quasi-condensate in the zero momentum state. 
Middle and bottom panels : As A is turned on, some particles 
leave the zero momentum state and go into higher momentum 
states. New peaks emerge at those momenta commensurate 
with the superlattice. These sharp peaks signal that the sys- 
tem is in a (gapless) superfluid state even if the superfluid 
density is very small (Fig. |HJ . 



the dimensional XY universality class when at in- 

teger densities the transition is driven by changing the 
on-site repulsion I n particular, in ID the 
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FIG. 10: The momentum distribution function for p = 0.5. 
Top panel : For A — 0, there is no superlattice potential and 
the distribution looks similar to Fig. [5] corresponding to a 
quasi-condensate. Middle and bottom panels : For A = 2 and 
A = 3, the density is commensurate with the superlattice 
leading to an insulating state. The momentum distribution 
very broad due to the absence of quasi-condensation. No ad- 
ditional peaks are observed at momenta conmensurate with 
the superlattice. These two features signal the presence of a 
(gapped) insulating state. 
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Before discussing the general case, t ^ 0, U ^ 0, and 
A 7^ 0, a useful step is to consider first the atomic limit 
t = 0. In this case the Hamiltonian is diagonal in the oc- 
cupation number basis, and the ground state is obtained 
by filling the system in a way which minimizes the en- 
ergy, depending on the competition between U and A. 
The on-site repulsion on one hand tends to avoid mul- 
tiple occupancies, leading to Mott insulators at integer 
densities. On the other hand, the modulating potential 
tends to trap the particles into its minima, leading to a 
density profile which reflects the modulation. In particu- 
lar, for T = 2, the modulating potential tends to impose 
a charge density wave (CDW) which results in alterna- 
tively highly and weakly occupied sites. 

For a given on-site repulsion U, we start to fill the 
system from p = by putting the particles alone on 
low energy sites until p = |. Each time a particle is 
added, the energy decreases by steps of —A, resulting 
in a chemical potential p = —A. Adding more particles 
leads to a competition between U and A. If A < U, the 
energy is minimized by putting the new particles on high 
energy sites, increasing the total energy by steps of A. 
On the contrary, if A > U, the on-site repulsion will not 
avoid double occupancies and the energy will increase by 
steps of 2U — A. As a result the value of the chemical 
potential from p = ^ to p = 1 will be p = min(A, 2U—A). 

Considering all possibilities of filling the system al- 
lows us to draw the phase diagram of the atomic limit 
in the (p/U, A/U) plane (Fig. I12f> . Regions labeled as 
"Mott" refer to configurations where the density p is 
integer with a uniform profile. The structure factor 




0.00 0.10 0.20 0.30 0.40 
t/U 



FIG. 11: (Color online) The phase diagram of soft-core 
bosons in the uniform case A = 0. 
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S(k) = {l/L 2 )J2jj> (njfij')e %k{ ~ ] J ^ then presents only 
one peak in k = with 5(0) = p 2 . Regions with label 
"CDWn" refer to staggered phases where the difference 
between the density on low energy sites and high energy 
sites is n. Phases CDWn with n even occur for integer 
densities, and phases with n odd for half-integer densi- 
ties. They have a signature in the structure factor which 
results in the development of a peak for k = 7r, with 
S(tt) = £. 



O Empty site 

— Particle on low 
energy site 

m Particle on high 
energy site 




A/U 



FIG. 12: (Color online) The phase diagram of soft-core 
bosons in the atomic limit t = (see text for details). 



B. General case 

For our QMC computations, we used the World Line 
algorithm. 24 ' 25 When turning on the modulating poten- 
tial, it is interesting to look first at the case with period 
T — 2. At finite, but large, on-site repulsion one can ex- 
pect to obtain incompressible regions for the same frac- 
tional densities as the hard-core case. This is because a 
small 1/U acts like a perturbation to the noninteracting 
spinlcss fermion Hamiltonian^ i.e., it should not change 
the nature of the phases present for infinite U. Moreover, 
for low values of A and sufficiently large U, we can expect 
incompressible regions with p — n (n being an integer), 
since this time it is A which acts as a perturbation on the 
translationally invariant boson Hubbard model (Eq. l(T|l). 
Figure El shows the density p as a function of the chem- 
ical potential p for T = 2, A = 2, and different values of 
U. The slopes of these curves, j^, are proportional to the 
isothermal compressibility. As a result, any discontinuity 
of (i (a gap) corresponds to a vanishing compressibility 
(presence of a plateau) . Starting with [7 = 1, we can see 
that such a "band" discontinuity occurs for p = |, as 
expected. As we increase U, the gap at p — h becomes 
larger and eventually a "Mott" gap opens at p = 1 as 
manifested by the plateau in Fig. ^] This corresponds 
to the first lobe in Fig. ^] Thus, the soft-core system 




FIG. 13: (Color online) The density of particles as a function 
of the chemical potential for T = 2, A = 2, and different 
values of the on-site repulsion U. 



is able to reproduce properties of the hard-core case (a 
gap at p = i), and properties of the uniform soft-core 
model (a gap at p = 1). In addition, a new gap for 
p = § appears starting from U — 3. Increasing the on- 
site repulsion further leads to the presence of a gap at 
p = 2, corresponding to the second lobe of Fig. EI Sim- 
ulations show that gaps appear also for other integer and 
half integer densities (p — §, §, |, • • •) with strong on-site 
repulsion. 

One can wonder how our QMC calculations are rele- 
vant to the zero temperature limit. Our algorithm de- 
termines the superfluid density by extrapolating to zero 
frequency the Fourier transform of the pseudo-current 
correlation function (j(r)j(0))m^ This gives the mean 
square value of the winding number, which is related to 
p s as defined by Eq. Ijl0|l £L The advantage of computing 
p s this way is that the measured value does not suffer 
from finite size effects (wc have considered L > 20), and 
can give the value of p s relevant to the zero temperature 
limit using an inverse temperature (3 not too large (we 
have taken (3 > 16). This is shown in Fig. 1141 which dis- 
plays exact analytical and numerical results of p s in the 
hard-core limit for O = and L = oo, and QMC compu- 
tations at finite temperature {(3 = 16), and finite system 
size (L — 20) as a function of the filling p. The data are 
in quite good agreement, and we can then expect this to 
hold in the soft-core case. 

In Fig. 1151 we show the superfluid density p s in the 
soft-core case as a function of the density p for A = 2, 
and T = 2. It is to be compared with Fig. 1131 and 1141 
Several things are apparent. Even at large U the super- 
fluid density is not anymore symmetric around p = 0.5 
showing the absence of the particle-hole symmetry of the 
hard-core case. The insulating phases (p s = 0), present 
at commensurate fillings for large U, start to disappear 
with decreasing U. Finally, for the smallest repulsive in- 
teraction we show in Fig.^]([/ =1), the only insulating 
phase occurs at p s = 0.5. The overall behavior is the one 
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FIG. 14: (Color online) Superfluid density p s as a function of 
density p. Comparison between exact analytical and numer- 
ical results at zero temperature in the thermodynamic limit, 
and an extrapolation using finite size and temperature QMC 
computations, for the uniform hard-core case, and for a case 
with A = 2 and T = 2. 




U/t 



FIG. 16: (Color online) The gap as a function of U for A = 
2, T = 2, p = i, 1, and several number of supercells. As 
the number of supercells increases, the "lattice" gap due to 
the finite size of the system decreases, showing evidence of a 
critical value of the on-site repulsion U for which the true gap 
vanishes. 




FIG. 15: (Color online) The superfluid density as a function 
of the density of particles, for T = 2, A = 2, and different 
values of the on-site repulsion U. 



one would expect after Fig. 1131 

As inferred from the results above, with decreasing on- 
site repulsion U one reaches a quantum critical point U c 
for which the gap at commensurate filling vanishes, lead- 
ing to a superfluid phase. This can be better seen in 
Fig. 1161 which shows the gap as a function of U for p = | 
and p = 1, and several system sizes. As U is lowered, the 
gap decreases but does not vanish completely if the sys- 
tem size is not large enough. This is due to finite size 
effects produced by the lattice gap. As the system size 
increases the lattice gap decreases and vanishes in the 
thermodynamic limit. Hence, in order to obtain a phase 
diagram one needs to do an extrapolation to check that 
the calculated gap does not depend on the system size. 

Runs similar to those presented in Fig. 1161 allowed us 
to obtain the phase diagram for soft-core bosons in the 



plane p/U vs t/U (Fig. I17|) . There one can see that the 
lobes at integer fillings are very similar to the ones of 
the homogeneous case depicted in Fig. ^] On the other 
hand, new lobes appear at n = 1/2,3/2,- • •. The first 
extends to rather low values of U for the case A = 2 
shown in Fig. 




_3 I i i i i i i i i i i i i i i i i i i i I 

0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 

t/U 

FIG. 17: (Color online) The phase diagram of soft-core 
bosons for A = 2 and T = 2 in the (p/U, t/U) plane. 

In order to establish a connection with the atomic 
limit, it is useful to draw a phase diagram in the 
(p/U, A/U) plane, as in Fig. 1121 for a fixed on-site re- 
pulsion U — 8t. This is done in Fig. 1181 where we use the 
notations "Mott" and "CDWn" to denote phases which 
have a resemblance with the ones in Fig. There are 
some similarities with the atomic limit. For example, the 
region with p = \ starts with p/U = for A/U = and 
has a maximum value of the chemical potential close to 
p/U = 1 for A/U = 1, then goes down and crosses the 
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axis p/U — at A/U « 2. In the region with p = 1, we 
have a transition from a Mott to a CDW2 at A/U = 1, 
corresponding to the critical point Pi. The same is true 
for p = 2 where we have a transition from a Mott to a 
CDW2 at A/U — 1, corresponding to the critical point 
P2. The critical point P3 is also present for p = |, and 
corresponds to a transition from a CDW1 to a CDW3. 
The main difference between the case U = 8 and the 
atomic limit t = is that the insulating regions with 
commensurated fillings are separated by incommensurate 
superfluid regions. In addition, even the Mott insulating 
phases exhibit a modulation in the density, in contrast to 
the constant density in the usual homogeneous case. This 
can be seen in the plots of the structure factor (Fig. IT§|> . 
which also signal clearly the transition between the dif- 
ferent phases shown in Fig. 1181 




0,25 0,5 0,75 1 1,25 1,5 1,75 2 2,25 2,5 



A/U 

FIG. 18: (Color online) The phase diagram of soft-core 
bosons for U = 8 and T = 2 in the (p/U, A/U) plane. 

As U and A are decreased no extrapolation is possible 
from the analytical results in the hard-core and atomic 
limits. QMC simulations are thus essential to understand 
this region. At commensurated fillings we find that the 
competition between U and A can drive the system su- 
perfluid over a finite range of values of U and A in be- 
tween Mott and CDW insulating phases. This can be 
seen in Fig. QUI where we have plotted the phase diagram 
for p = 1 and T = 2. Our results for intermediate val- 
ues of U and A not only contrast with the atomic limit 
case where no intermediate phase is present, but also with 
studies of a similar model for fermionic systems^ We are 
refering to the fermionic Hubbard model with an addi- 
tional T = 2 potential, also known as the Ionic Hubbard 
model. In this model an intermediate phase was also ob- 
served between the Mott insulating and band insulating 
phases. However, in the fermionic case the intermediate 
phase turned out to have a finite one particle gap^& while 
we find it to be gapless (superfluid) in our soft-core boson 
case. 

There are interesting qualitative and even quantitative 
analogies between the T — 2 phase diagram considered 
here, in which a Mott phase competes with a CDW phase 
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FIG. 19: (Color online) The structure factor S(k) shows a 
peak for k = n when the density profile is staggered. The 
value S(n) allows to detect the height of the density steps. In 
the atomic limit, the previously defined CDWn phases corre- 
spond to S(n) = Here, in the general case, these phases 
are recovered when the ratio A/U is away from the critical 
points Pi , P2 , P3 , Pi ■ 



driven by the one-body superlattice potential A and the 
phase diagram of the extended Hubbard model where 
CDW correlations arise from a near neighbor interaction 
In both cases, a superfluid region extends along a 
strong coupling line out to U f» 6t, and the superfluid 
extends to arbitrarily large V or A at U — > 0. 




0123456789 



U/t 

FIG. 20: (Color online) The phase diagram of soft-core 
bosons for p — 1 and T = 2 in the (A/t,U/t) plane. 



A behavior similar to the T — 2 case can be expected 
to hold for larger periods of the superlattice. In the ap- 
pendix we present some results for the case T = 6 that 
support this conclusion. 
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IV. THEORETICAL ANALYSIS OF T = 2 
MODEL 

In the atomic limit (t — ► 0), the sites are decoupled and 
described by atomic Hamiltonians Hj = (Vj—fi) h + Un 2 , 
where Vj = A cos(27rj/T). As (i, A, or U is var- 
ied, there are a series of first order transitions between 
number eigenstates, occurring when either the particle 
(A p = E(n + 1) - E{n)) or hole (A h = E(n - 1) - E{n)) 
gaps collapse, where 

A p = 17(1 + 2n) + Vj - ii (14) 
A h = U(l - 2n) -Vj+fi . (15) 

Thus, the atomic state | K) is the local ground state for 

U(2n - 1) < n - Vj < U(2n + 1) . (16) 

Note that A p (n - 1) = -A h (n). 

Consider now the T = 2 case, for which Vij = +A and 
Vij+i = —A. In the atomic limit, the eigenstates are of 
the form 

i^>*°>=n%^%pH o > ( i7 ) 

^ V e • V o • 

The state | K c , K Q ) is the ground state provided all four 
gaps A®, A°, A|, and AJJ are positive, where 

Ap\ h = U ± {2UK C + A — n) (18) 
A°_ h = *7 ± (2C/X — A — n) . (19) 

These four inequalities define the colored rectangular re- 
gions in Fig. El (see also the detail in Fig. 1211 1 




\K.-l,K-l) <!■> 



FIG. 21: (Color online) Detail of phase diagram of the T = 2 
model in the t — limit. The gaps are defined relative to the 
shaded region. 



Let us now investigate the case where t/A is very small. 
We begin with the atomic state \K e ,K ). We first as- 
sume that A° and Ag are always positive and much larger 
than A£ and A°, which puts us in the right corner of the 
shaded region in Fig. Focusing on the lowest-lying ex- 
citations, which are holes (particles) on even (odd) sites, 
we arrive at the effective Hamiltonian 

H = -t Y,( h hph+i + h hph-i+ R - c -) ( 2 °) 

3 

+ e ( A h h h k + a p ph+i P23+1) 

+ U Y1 ( h h h h k 2j h 2j + Plj+l Ptj+l P2j+1 Pij+l) ' 
3 

where h\ - creates a hole (destroys a boson) on site 2j and 

p\j+i crea tes a particle (creates a boson) on site 2j + 1. 

The hopping integral is now t w ty/K c {K a + 1). In the 
continuum limit, the coherent state Lagrangian density 
is 

C = h(d T + A h )h + p(d T + A p )p-2b- 1 t(hp + hp) 

+ i bt (d x hd x p + d x h d x p) 

+U{hhf + U{ppf , (21) 

where b is the unit cell length, and where we abbreviate 
A h = A° and A p = A°. 

Suppose A^ ^ A°, and further assume that U <C t is 
very weak. Then we can integrate out the hole states in 
the low energy limit, obtaining the effective action 

f, f fdk ( 4? Pb 2 k 2 \ n . . 

s = rxlTA 1 -^^)^ (22) 

/ At 2 t 2 b 2 k 2 \ -r . . 2 1 

where U/U = l+0(t 2 / A 2 ). This predicts a transition to 
a compressible phase when the effective p gap collapses, 
which occurs at 

A p A h = P . (23) 

Alternatively, we can write a trial ground state for the 
model, where h(x)\ |*) = h\^) and p(x)\ |*) = 
The variational energy is minimized when the product 
p h is real, and we may then assume both p and h are 
real. Varying with respect to the amplitudes p and h, 
one obtains the coupled nonlinear equations 

A p p + 2Up 3 -2th = (24) 
A h h + 2Uh 3 - 2tp = . (25) 

For A p > and Ah > and A p Ah > t 2 , the solution is 
p = h = (i.e. a CDW state). For A p < and A h < 
and A p Ah > t 2 , the solution, if we assume t is weak, is 
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p ~ y/—A p /2U and h ~ — \/— A^/2U, which is the best 
the coherent state description can do in approximating 
the neighboring CDW state \K e -l, K Q + 1) . If A p A h < 
i 2 in cither of these cases, the state is compressible. 

If U S> Apj^i, we may assume that the number of 
holes or particles on each site is either zero or one, and 
map the problem back onto the soluble U — oo model. 
We stress that U does enter into the formulae for Ap' 1 ^, 
so this is not entirely trivial. Again we focus on the right 
corner of the shaded region in Fig. |^ The \K e ,K ) 
state is represented as | Tl), i.e. a state where the 'spin' 
on each even (odd) site is polarized up (down). The 
effective S — \ Hamiltonian is then 

H = —t (s£ S n+1 + S n S^ + ij 

n 

+A h ^Sf i -A p £Sf J . +1 . (26) 

3 3 

Solving again via Jordan- Wigner fermionization, one has 
two energy bands in the reduced zone fee [ — -|, ^1 with 
dispersions 

E ± (k) = i(A h - A p ) ± i^/(A h + A p )2 + At 2 . (27) 

This leads to the following phase classification: 

A h > , A p > : incompressible ~ | tl) (28) 

A h < , A p < : incompressible ~ | J.f) (29) 

A h < , A p > : |TT) if |A P A h | > t 2 (30) 

A h >0 , A p <0 : Hi) if |ApA h | >t 2 . (31) 

Thus, the incompressible phases occur just above and 
below the nodes of the phase diagram, where four atomic 
phases meet. This is in qualitative agreement with the 
numerically obtained phase diagram in Fig. 1121 

V. CONCLUSIONS 

The interplay between particle-particle interactions 
and a one body potential on the phases of correlated 
quantum systems is a fascinating and complicated ques- 
tion. In the case of a random one-body potential, the 
key issue is whether interactions between electrons can 
cause an Anderson insulator to become metallic, par- 
ticularly in two dimensions, a question which has not 
been definitively resolved either experimentally or theo- 
eretically. The effect of interactions on a band insulator 
has recently been explored for the ionic Hubbard model 
in one dimension, with the interesting suggestion that, 
as in the extended Hubbard model where CDW corre- 
lations arise from interactions, there is a bond-ordered 
wave phase in a region where spin and charge order cor- 
relations are in a delicate balance. 

This paper has provided a careful examination of the 
effect of correlations on boson systems in a superlattice 



potential in one dimension. This is a particularly inter- 
esting case to explore, since the hard-core limit connects 
to the fermion problem. Indeed, as we have shown, the 
band insulating behavior present in the hard-core case 
seems to persist when U is finite, even though the bosons 
can now multiply-occupy the sites, and one no longer has 
concepts like the Pauli principle, a Fermi-surface, etc... 
which are key ingredients to the usual picture of a band 
insulator. Furthermore, we have shown that at integer 
fillings the transitions between insulating phases, driven 
by changing the ratio on-site repulsion / strength of the 
superlattice potential, can produce intermediate super- 
fluid regions. This is something that could be tested in 
experiments with ultracold gases on optical lattices since 
in these intermediate regions coherence is enhanced in 
contrast to the insulating Mott and CDW phases. 
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APPENDIX: A: LOCATOR EXPANSION FOR 
THE GREEN'S FUNCTION 

In the U = oo case, several analytical results for 
the energy width of the T bands may be obtained us- 
ing a locator expansion for the Green's function. For 
each n S {0, . . . , T — 1}, construct the quasi-momentum 
eigenstate \n{k)) = (T/L) 1 / 2 J2j e ik ^ n + j7 "> \n + jT). In 
this basis, the Hamiltonian matrix is 

fit* = Acos^Cx-ie-^e.n-r 

-te^S T n , n+1 . (A.l) 

Consider a ring of T sites, which represent the con- 
secutive states \n(k)). We now compute the Green's 
function G{E) = (E - H)" 1 in this basis, using a lo- 
cator expansion in powers of the hopping t. The bare 
(t = 0) Green's function is diagonal in this basis, with 
G° m = [E — Ell)^ 1 . Consider first a state which is non- 
degenerate in the atomic limit, i.e. n — for T odd, 
and n = and n = T/2 for T even. Let be the 

self-energy contribution from all paths which start and 
end at n but do not contain n as an intermediate state, 
and which have zero net winding number (E°) or wind 
once clockwise (£~) or wind once counterclockwise (S + ). 
Summing the perturbation series for G nn yields 

Q-l _ irfi \-l _ yO _ y+ _ y- ( A n\ 

^ nn V nn/ nn nn nn ' \ / 

The full i-dependence of the energy level E n (k) is now 
given by the solution to the equation 

E = E° n (k) + E°„ (E) + £+„ (E) + E" n (E) , (A.3) 
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which results in a pole in G nn (E). Since the self energy 
terms do not contain any factors of G°„, we can to a first 
approximation set E = E® therein. The self energy S°„ 
contains infinite orders in t, but it does not accrue any 
phase, since it corresponds to paths of zero net winding. 
Therefore, it leads to a fc-independent shift of the atomic 
energy levels. The fc-dependence enters through the self 
energies S^„. Since we are interested in the small t limit, 
we evaluate the lowest order contributions: 



KuM) = ^ (fe+c)T 't T n 



i 



EO 



(A.4) 



with £-„=(£+„)*. Thus, 



B„ t 



E n (k) = E° n + AE n (t, A) - -^-j- cos(Tfc + T() + . . . , 

(A.5) 

where B n is a constant. We now read off that the band- 
width is of order tF jA T ~ x . 

Next consider the case of degenerate atomic levels n 
and n = T — n. We define the self energies E°„ and 
£5 fi as the self energy contribution from all paths start- 
ing and ending at n or h and which contain neither n 
nor n as intermediate states. We additionally define T^n 
as the self energy contribution from all paths starting at 
n and ending at n, circulating clockwise (£~) or coun- 
terclockwise (£ + ), and which contain neither n nor n as 
an intermediate state. A corresponding definition holds 
for £^„, from which it follows that £^„ = (E^)*. The 
locator expansion for G nn = Gnn may be summed: 



sy-l _ (pQ \-l 
^nn V^nn) 



nn 



y—_ 

nn 



(A.6) 



Thus, the degenerate levels split, and are given by solu- 
tions to the equations 

E = E a n {k) + £°„(£0 ± | E+ n (E) + E~ n (E) \ . (A.7) 

Again, is fc-independent, and evaluating the k- 

dependent self energies S^ fi to lowest order in t, we ob- 
tain (with 1 < n < 



^nn( E n) 



te i(k+0 



, T-2r, 



nj=„+r (k-e?) 



(A. 



(_ te -*(fc+C))< 



t 



T-2r, 



^T-2n-l 



3 -i(T-2n)(fe+C) 



j.2n 

n p »(2n)(fe+C) 



(A.9) 



if 1 < n < -j, and 



^tn(En) + ^nn 



C n t? 



(El) 

2n Y) 



(A.ll) 



2 ii 



,((T-4n)(fc + 0) 



j^T-2n-l ' J{2n~l 

j < n < -j. Thus the bandwidth T n as well as the 
inverse effective mass (m*) _1 for the bands arising from 
the atomic levels I n ) and I T — n ) scales as 



if 



/ l \ max (2ri-l,T-2n-l) 



(A.12) 



where 1 < n < -j. As we shall see, this power law be- 
havior also governs the scaling of the superfluid density 
with t/A in the large A limit. 

Further consideration of Eq. (|A.7|l shows that the de- 
generate zero energy states at n = p and n = 3p for 
T = Ap do not shift for k = £ — 0, a statement valid 
to all orders in t/A. The reason is that the self energy 
contributions S° p and T,pp + T,~p each vanish at E = 0, as 
a consequence of the fact that E^ + . = — Jg£ ., and hence 
for every path contributing to these two self energy con- 
tributions there exists a path with equal and opposite 
amplitude, resulting in a cancellation in the locator ex- 
pansion. 



APPENDIX: B: T = 6 CASE FOR SOFT-CORE 
BOSONS 

Figure [221 displays the density p as a function of the 
chemical potential p, for A — 2 and U = 8. As for the 
hard-core case, gaps appear for fractional fillings p — 
i, |, |, |, |. The integer fillings p = 1, 2, • • • are also 
insulating, as expected from the uniform soft-core case. 
Again, extra gaps also appear for fractional fillings with 
p>l. 

The inset in Figure [221 shows the density profile for 
[7 = 4 and {7=8. As for the T = 2 case, the local 



Thus, 




^tn( E n) + ^nn( E l) \ 

D n t 2n . C„i T " 2 " 



(A.10) 



j^2n-l 



^T-2n-l 



,((T-4n)(fc + C)) 



FIG. 22: (Color online) The density of particles as a function 
of the chemical potential for T — 6, A = 2, and U — 8. 
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density does not stick to p L = 1. Finally, Fig. 1231 shows 
the superfluid density which vanishes for almost each fill- 
ing commensurate with the superlattice. There are some 
exceptions for p — | , |, | due to the small values of the 
gap (see Fig. 0^1 . 




1/6 1/3 1/2 2/3 5/6 1 7/6 4/3 3/2 5/3 11/6 2 13/6 7/3 5/2 
P 

FIG. 23: The superfluid density as a function of the density 
of particles, for T = 6, A = 2 and U = 8. 
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